!------------------------------------------------------------------------------
! TJU/Department of Mechanics, Fluid Mechanics, Code START
!------------------------------------------------------------------------------
!> @Author: Jianxin Liu shookware@tju.edu.cn
!> @Date: 2024-08-13 14:27:01
!> @LastEditTime: 2024-08-13 14:27:03
!> @LastEditors: Jianxin Liu shookware@tju.edu.cn
!> @Description:
!> @FilePath: /START_NS_3D/src/finite_difference_scheme/submodule_finite_difference_scheme.f90
!> @Copyright (c) 2024 by Jianxin Liu email: shookware@tju.edu.cn, All Rights Reserved.
!------------------------------------------------------------------------------
submodule(module_finite_difference_scheme) submod_fdm_4center

contains

   module subroutine CreateFDMScheme_4Center(this, length, stencil_width, flg_bounds)

      implicit none
      class(type_fdm_scheme_4center), intent(inout) :: this
      integer, intent(in) :: length
      integer, intent(in) :: stencil_width
      logical, intent(in) :: flg_bounds(2)

      integer :: i
      integer :: sw, iend

      call this%set_parameters(length, 2, flg_bounds)

      call this%D%create(this%length)
      call this%D%setup()
      call this%D%diagm(-2, [(1.0d0/12.0d0, i=0, this%length - 3)])
      call this%D%diagm(-1, [(-8.0d0/12.0d0, i=0, this%length - 2)])
      call this%D%diagm(1, [(8.0d0/12.0d0, i=0, this%length - 2)])
      call this%D%diagm(2, [(-1.0d0/12.0d0, i=0, this%length - 3)])
      sw = this%stencil_width; iend = this%length - 1

      if (this%flg_bounds(1)) then
         call this%D%set_value(sw + 0, [(sw + i, i=0, 2)], [-3.0d0/2.0d0, 2.0d0, -1.0d0/2.0d0])
         call this%D%set_value(sw + 1, [(sw + i, i=0, 3)], [-2.0d0/6.0d0, -3.0d0/6.0d0, 6.0d0/6.0d0, -1.0d0/6.0d0])
      else
         call this%D%set_value(0, [(i, i=0, 2)], [-3.0d0/2.0d0, 2.0d0, -1.0d0/2.0d0])
         call this%D%set_value(1, [(i, i=0, 3)], [-2.0d0/6.0d0, -3.0d0/6.0d0, 6.0d0/6.0d0, -1.0d0/6.0d0])
      end if
      if (this%flg_bounds(2)) then
         call this%D%set_value(iend - (sw + 1), [(iend - (sw + 1 + i), i=-1, 2)], &
         & -[-2.0d0/6.0d0, -3.0d0/6.0d0, 6.0d0/6.0d0, -1.0d0/6.0d0])
         call this%D%set_value(iend - (sw + 0), [(iend - (sw + i), i=0, 2)], -[-3.0d0/2.0d0, 2.0d0, -1.0d0/2.0d0])
      else
         call this%D%set_value(iend - 1, [(iend - (1 + i), i=-1, 2)], &
         & -[-2.0d0/6.0d0, -3.0d0/6.0d0, 6.0d0/6.0d0, -1.0d0/6.0d0])
         call this%D%set_value(iend - 0, [(iend - i, i=0, 2)], -[-3.0d0/2.0d0, 2.0d0, -1.0d0/2.0d0])
      end if
      call this%D%assembly()
      if (this%flg_bounds(1)) call this%D%zero_rows_columns([(i, i=0, sw - 1)], 0.0d0)
      if (this%flg_bounds(2)) call this%D%zero_rows_columns([(iend - i, i=0, sw - 1)], 0.0d0)

   end subroutine CreateFDMScheme_4Center

   module subroutine DiffFDMScheme_4Center(this, f, df)
      implicit none
      class(type_fdm_scheme_4center), intent(inout) :: this
      class(type_vector_base), intent(in) :: f
      class(type_vector_base), intent(inout) :: df

      integer :: ierr
      call this%D%matmult(f, df)

   end subroutine DiffFDMScheme_4Center

   module subroutine DestroyFDMScheme_4Center(this)

      implicit none
      class(type_fdm_scheme_4center), intent(inout) :: this

      call this%D%destroy()
      this%length = -1
      this%stencil_width = -1
      this%flg_bounds(2) = .FALSE.

   end subroutine DestroyFDMScheme_4center

end submodule submod_fdm_4center
